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Carbon sequestration in geologic reservoirs is an important approach for mitigating greenhouse gases 
emissions to the atmosphere. This study first develops an integrated Monte Carlo method for simulating 
CO2 and brine leakage from carbon sequestration and subsequent geochemical interactions in shallow 
aquifers. Then, we estimate probability distributions of five risk proxies related to the likelihood and volume 
of changes in pH, total dissolved solids, and trace concentrations of lead, arsenic, and cadmium for two 
possible consequence thresholds. The results indicate that shallow groundwater resources may degrade 
locally around leakage points by reduced pH and increased total dissolved solids (TDS). The volumes of pH 
and TDS plumes are most sensitive to aquifer porosity, permeability, and CO2 and brine leakage rates. The 
estimated plume size of pH change is the largest, while that of cadmium is the smallest among the risk 
proxies. Plume volume distributions of arsenic and lead are similar to those of TDS. The scientific results 
from this study provide substantial insight for understanding risks of deep fluids leaking into shallow 
aquifers, determining the area of review, and designing monitoring networks at carbon sequestration sites. 

Geologic carbon sequestration (GCS) in deep reservoirs is considered to be a viable long-term solution for 
mitigating greenhouse gases in the atmosphere\ But there are substantial regulatory and technical 
challenges that need to be overcome before GCS can be widely deployed. A major regulatory concern 
that has garnered much research attention is the potential leakage of CO2 and brine from deep storage reservoirs 
to shallow groundwater resources through preferential pathways, such as faults, wells, and local high-permeab- 
ility zones in caprocks. The leaked CO2 and brine may cause a decrease in groundwater pH, an increase in total 
dissolved solids (TDS), and a potential mobilization of trace metals from aquifer materials into groundwater^"^". 
While the potential risks to shallow aquifers at GCS operations are known, no defensible methodology for 
quantifying the risks has been reported to date^'^. 

As part of the National Risk Assessment Partnership (NRAP) project^\ we are developing approaches that can 
be used to quantify long-term risks at GCS sites. Current efforts are focused on analysis of risk proxies (i.e., 
impacts and their likelihoods) that will ultimately be combined with consequence analyses to create rigorous risk 
assessments. NRAP utilizes an Integrated Assessment Model (lAM) to predict long-term performance of GCS 
sites. The lAM is developed using a system modeling approach that captures the physical and geochemical 
processes that occur within and among various sub-components of the sequestration sites, e.g., from sequest- 
ration reservoir to leaky well, then to underground sources of drinking water (USDW) and/or the atmosphere. As 
an integrated model, the I AM simulates the coupling processes from the point of CO2 injection to the interactions 
among CO2, brine, groundwater, and aquifer materials in case of a leak. Coupling the complex sub-component 
models involves significant computational costs, since the inherent heterogeneities and uncertainties in geologic 
systems require modeling a large number of realizations to bracket and quantify uncertainties associated with 
likelihoods, impacts, and risks. To meet the requirements of modeling coupled processes and performing large 
number of realizations, NRAP is developing and using reduced order models which are computationally efficient 
versions of the full process models that capture the underlying physical and chemical processes at reduced 
computational costs^\ 

This study demonstrates the developed methodology by using the highly heterogeneous unconfined Edwards 
aquifer as an example (Texas, USA). This aquifer was used because there is considerable data available, and 
because carbonate aquifers represent an important class of drinking water sources in the USA. Furthermore, the 
Edwards aquifer overlies potential sequestration sites including oil and gas reservoirs^^ The Edwards aquifer is 
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P6 


Max. CO2 pressure (MPa) 1 0 


40 


P7 


Max. CO2 saturation 0.1 


1.0 


P8 


Wellbore permeability (m^) 1 0"^^ 


10-^2 


P9 


CI concentration (mol/l) 0.5 


5.4 



one of the most productive Karst aquifers in the world. It serves the 
diverse agricultural, industrial, recreational, and domestic needs of 
almost two million users in south central Texas^^. The total ground- 
water withdrawal from the Edwards aquifer in 2000 was 2.8 X 
10^ mVd of which withdrawal for public supply was about 1.6 X 
10^ m^/d. The Edwards aquifer consists of an underground layer of 
porous, fractured, honeycombed, water-bearing limestone rock with 
highly heterogeneous structures. In the recharge zone, the aquifer is 
unconfined. In the transition and discharge zones, it is an artesian 
aquifer confined by a layer of very low permeability sediment. 

We select an unconfined zone of the Edwards aquifer as an 
example in order to evaluate the impact of the leaked CO2 and brine 
on the USDW and the potential CO2 discharge rate into the atmo- 
sphere (Note that other NRAP groups use the confined High Plain 
aquifer as an example to cover different types of shallow aquifers in 
the USA^^). The study area is located north of San Antonio, where the 
unconfined carbonate aquifer has a total thickness between 100 and 
200 m, with a mean thickness of 150 m. The groundwater flows from 
northwest to southeast with a hydraulic gradient around 0.00087 and 
the permeability structure of the unconfined carbonate is highly 
heterogeneous^'*. By reviewing the previous study of the permeability 
data^^"^^ we summarize the statistical and spatial correlation para- 
meters such as mean, variance, and integral scale in Table 1. The 
vertical integral scale is assumed to be linealy correlated to the 
horizontal integral scale with a factor of 0.01. The existing permeab- 
ility data from Lindgren et al.^^"^^ are used as pilot-point data for 
geostatistical simulations. A pilot-point-based Gaussian simulation 
method is modified from GSLIB^^ to a geotatistical model GEOST^^'^* 
for simulating permeability fields. PSUADE^^ and GEOST are 
coupled to generate nearly 500 realizations of the process model with 
Latin Hypercube sampling and geostatistical modeling for an 



integrated Monte Carlo simulation of CO2 reactive transport in the 
unconfined carbonate aquifer. 

For each realization the uncertain parameters are sampled first 
and then the heterogeneous permeability fields for the unconfined 
aquifer are simulated with the pilot-point-based Gaussian method. 
Figure 1 shows one realization of the generated permeability fields in 
three dimensions. The model y axis is set along the groundwater flow 
direction and the lateral flow rate is varied with the sampled per- 
meability and porosity. The numerical model size is 5000 X 8000 X 
150 m with 64315 elements. The possible CO2 and brine leakage 
point is set at the upstream with a coordinate of (2500, 7000, 0), 
where the grid is highly refined with minimum grid sizes Ax, Ay, 
Az of 2, 4, and 1 m, respectively. Away from the leakage point, the 
numerical grids become coarse (Figure 1). 

The maximum CO2 pressure and saturation in the sequestration 
reservoirs and the wellbore permeability in the caprocks listed in 
Table 1 are three parameters we used to estimate the CO2 and brine 
leakage rates from the deep reservoirs to the shallow aquifer with a 
reduced order model developed as part of NRAP^^. The reduced 
order model computes CO2 and brine leakage source terms using 
CO2 pressure and saturation in sequestration reservoir and the well- 
bore permeability. One example of the CO2 and brine leakage func- 
tions is shown in Figure 2, in which the CO2 injection time and rate 
in the deep reservoirs were assumed to be 50 years and 5 million 
tons per year, respectively. When CO2 and brine leaks into the 
unconfined carbonate aquifer, there could be very complex CO2- 
induced geochemical reactions in the aquifer, including aqueous 
equilibrium reactions (or acid-base and complexation reactions), 
trace metal adsorption or desorption, and mineral dissolution or 
precipitation* *'*^'^^ After reviewing all those reactions discussed by 
Bacon^^ and Zheng et al.*^ for carbonate aquifers such as Edwards 
aquifer, we incorporated the reactions, which have the most im- 
pact on the pH plume development, into our reactive transport 
simulations: 

(1) CO2 gas dissolution in water to reduce pH, 

C02(gas) ^C02(aqueous) + HlO^HCO^ + H + , 

(2) Aqueous equilibrium reactions and the equilibrium constants 
at 25°C, 

HCO^+H+ =H2C03 logK= -6.3, 
CO^"+2H+ =H2C03 logK=-16.6, 




Figure 1 | The simulated heterogeneous permeability (logm^) field in a 3-D view. 
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Figure 2 | The simulated CO2 and brine leakage rates in a leaky well from 
one realization. 

OH-+H+ =H20 logK= -14, 

(3) Calcite kinetic dissolution or precipitation, 

CaC03(soiid) + 2H+ = H2CO3 + Ca'+ . 



Calcite is the most plentiful mineral in the carbonate aquifer and 
can significantly buffer the pH. The forward rate for the kinetically- 
contr oiled calcite dissolution or precipitation at 25°C is 1.5 X 
10"^ mol/mVs^\ The specific surface area for calcite is 0.1 cm^/g. 
Dolomite is another common mineral in carbonate aquifers. In this 
study area the volume fraction of dolomite is very small relatively to 
that of calcite. We ignored dolomite in our geochemical calculations. 
We assumed that the brine is mainly made up of NaCl and treat its 
concentration as an uncertain variable, ranging from 0.5 to 5.4 mol/ 
L. The geochemical analysis conducted by Bacon^^ and Last et al.^^ 
indicates that As, Cd, and Pb concentrations in the brine are pro- 
portional to the NaCl concentration. The proportional ratios are 3.16 
X 10"', 3.16 X 10"' and 1.0 X 10"' on a molar basis for As, Cd, and 
Pb, respectively. We simulated the species CI, Na, and trace metals 
(As, Cd, and Pb) from the leaked brine as conservative tracers by 
ignoring the sorption reactions of the trace metals. Since sorption 
processes generally reduce the concentrations of the trace metals, our 
approach represents the conservative or worst-case scenarios for the 
trace metals. 

By examining the existing water quality data in the unconfined 
Edwards aquifer we establish baseline data sets and statistical pro- 
tocols for determining statistically significant changes between 
background concentrations and predicted concentrations (after 
CO2 and brine leakage) The background (or initial) concentra- 
tions of the major chemical components, as well as their EPA max- 
imum contaminant levels (MCLs) and no-impact thresholds, are 
listed in Table 2^^'^^. Note that the no -impact thresholds are deter- 
mined for pH, TDS, As, Cd, and Pb based on their standard devia- 
tions of the measured water quality data and are used to identify 
potential areas (or volumes) of contamination predicted by aquifer 
numerical models. This study quantitatively evaluates the following 
five volume-based risk proxies: 1) plume volume of pH less than the 
MCL or a no-impact threshold in the aquifer; 2) plume volumes of 
TDS larger than the MCL or a no-impact threshold; and 3) plume 
volumes of the three trace metals (e.g. As, Cd, and Pb) over their 
MCLs and no-impact thresholds. Additionally, we calculate the CO2 
discharge rate into the atmosphere (It should be noted that cur- 
rently there is no consequence threshold used for CO2 discharge 
rate to the atmosphere). 



Table 2 | Initial values, no-impact thresholds, and MCLs for each 
component 



Component 


Initial values 


No-impact thresholds 


MCLs 


Units 


Arsenic 


0.31 


0.55 


10 




Cadmium 


0 


0.04 


5 




Lead 


0.06 


0.15 


15 


^g/L 


pH 


6.9 


6.6 


6.5 


standard 


TDS 


330 


420 


500 


mg/L 



Results 

The integrated MC simulation has been developed in this study for 
assessing the plume volumes of pH, TDS, and trace metals (e.g. As, 
Cd, and Pb), and CO2 discharge rates into the atmosphere under a 
probability framework. In most of the MC realizations the shallow 
groundwater resources are degraded locally around leakage points by 
reduced pH, increased total dissolved solids, and trace metals mobi- 
lized from the injection reservoir. The simulated pH plumes are the 
largest of the risk proxies and extend from 100 to 1500 m along the 
flow direction. Perpendicularly to the flow direction pH plumes 
extend less than 500 m. Subsequent to the leakage free phase of 
CO2 gas is formed, which preferentially migrates upwards due to 
buoyancy and leads to formation of pH plumes to the water table. 
The amount of CO2 leaving the top of the aquifer ranges between 0 - 
98% of the CO2 leaking into the aquifer. 

The global sensitivity results indicate that the volumes of pH less 
than the MCL and the no -impact threshold are most sensitive to 
aquifer porosity and CO2 leakage rates, as well as the cumulative 
CO2 mass and aquifer mean permeability. The volumes of TDS larger 
than the MCL and the no -impact threshold are most sensitive to the 
cumulative brine mass, brine leakage rate, and aquifer porosity. The 
CO2 discharge rates (into the atmosphere) are positively correlated to 
the leakage rate from deep reservoirs. 

For all of the output variables pH plume volumes (measured with 
MCL or no -impact thresholds) have the largest mean and median 
while the trace metal Cd plume volumes are the smallest. When 
measured with no -impact thresholds, the plume volume distribu- 
tions of As and Pb are similar to those of TDS. Since the pH plume 
sizes are the largest among the five output variables, we will take it as 
one of the major risk proxies and use it as the reference to determine 
the area of review and to design the monitoring network in the 
unconfined shallow aquifer. 

Discussion 

The numerical results from the multi-phase reactive transport simu- 
lator FEHM^^ show that when CO2 leaks into the aquifer, there is a 
small fraction of CO2 gas (less 10%) dissolved into water, which leads 
immediately to a pH reduction at the leak points. The groundwater 
lateral flow transports the gaseous and dissolved liquid CO2 down- 
stream to extend the pH plumes. At the same time, the reduced pH 
enhances calcite dissolution to consume the increasing H^, which 
buffers the pH. Because the leaked gaseous CO2 has a much smaller 
density than water in the aquifer, the gaseous CO2 also transports 
upwards by buoyancy, which leads to pH plumes to extend vertically 
until the water table (see Figures 3 A, 3B and 3C). Finally the gaseous 
CO2 leaves the top of the unconfined aquifer into the atmosphere. In 
most of the MC realizations, the simulated pH < 6.5 plumes extend 
from 100 to 1500 m along the flow direction. Perpendicularly to the 
flow direction pH plumes extend less than 500 m (Figure 3D). 

Global sensitivity analysis. By using the Monte Carlo simulation 
results as the input for PSUADE^^, we conduct a global sensitivity 
analysis with the Main Effect method. At this step the three 
parameters that define the CO2 and brine leakage functions (e.g. 
maximum CO2 pressure and saturation in the squestration 
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Figure 3 | Simulated pH plumes in the aquifer at different times from a YZ plane view at 40 (A), 100 (B), and 200 years (C), and a XY plane view (D) at 
four different times. 



reservoirs and the wellbore permeability in the caprock) are 
converted to time- dependent leakage rates and the accumulated 
CO2 mass computed from the leakage functions for each 
realization. The results plotted in Figure 4 show that the plume 
volumes measured with MCLs and no -impact thresholds have the 
similar sensitivity to the input parameters, but different risk proxies 
(e.g. pH and TDS volumes) are sensitive to different parameters. For 
example, the volumes of pH less than the MCL and no-impact 
threshold are most sensitive to aquifer porosity and CO2 leakage 
rates, as well as the cumulative CO2 mass and aquifer mean 
permeability (Figure 4A). Since the aquifer porosity is linearly 
correlated to the volumes of pH plumes, it is most sensitive to that 
parameter. The volumes of TDS larger than the MCL and no-impact 
threshold are most sensitive to the cumulative brine mass, brine 
leakage rate and aquifer porosity (Figure 4B). The sensitivities of 
the volumes of the tracer metals (As, Cd, and Pb) less than their 



MCLs and thresholds are similar to those of the volumes of TDS. 
The CO2 discharge rate (into the atmosphere) is positively correlated 
to the leakage rate from the deep reservoirs. 

Probability measures for the risk proxies. By using the post- 
processing results of the original 460 (200 year) MC simulations, 
we conduct a statistical analysis of the output variables with their 
MCLs and no-impact thresholds. The 5^ 50'^ (median), and 95*^ 
percentiles of the outputs are computed at 20 different time steps. 
The statistical analysis results of the CO2 discharge rate to the 
atmosphere are shown in Figure 5. 

Since the MCLs and no -impact thresholds for pH and TDS are 
very close to each other (Table 2), the computed pH and TDS plume 
volumes with thresholds are slightly higher than their volumes with 
MCLs (Figures 5 A and 5B). The MCLs and no-impact thresholds for 
trace metals (As, Cd, and Pb) are quite different and their percentile 



(A) 



100 
80 
60 
40 
20 
0 



■ pHvolume-MCLs 

■ pH volume-Threshold 



ijjLii 



J- > .0*^ «< 



M?* (O" 




Figure 4 | Global sensitivity of the volumes of pH (A) and TDS plume volumes (B) to the 10 input parameters. 
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Figure 5 | The statistical analysis results for the 5 groundwater risk proxies and atmospheric leakage of CO2 with their 5th, 50th (median), and 95th 
percentiles. 



curves are also quite different (Figures 5C, 5D, and 5E). With high 
standard MCLs the plume volumes of the three trace metals are 0 at 
the 5*^ percentiles. The median plume volume of Cd is also 0 when 
measuring with its MCL (Figure 5D). The trace metal plume volumes 
measured with no-impact thresholds are significantly larger than 
those measured with MCLs. Surprisingly, the computed 5*, 50* 
(median), and 95* percentiles of the trace metal plume volumes with 
no-impact thresholds are very similar to the TDS volumes. Note that 
all these plume volumes still keep increasing slowly even though the 
CO2 and brine leakage decreases after 50 years. The time- dependent 
trend of the of the CO2 discharge rate to the atmosphere is similar to 
the CO2 leakage rate to the aquifer, but the discharge rate has a larger 
tail at the late time, which means that when the leakage rate (into 



aquifer) goes down after 50 years the gaseous CO2 in the unconfined 
aquifer still keeps releasing into the atmosphere with a gradually 
reduced rate. The 5* percentile of the CO2 discharge rate is 0 (a value 
of 10"^ in the log scale to represent 0 in Figure 5F). 

Figure 6 shows histograms of the plume volume distributions of 
the 5 risk proxies at 200 years. Note that the log volume of 0 repre- 
sents no plume. Among all these 5 risk proxies pH plume volumes 
(both with MCLs and thresholds) have the largest mean and median 
while the Cd plume volumes are the smallest. For As and Pb, 33% of 
the 460 realizations do not exceed the MCL threshold and 80% of the 
realizations do not exceed the MCL for Cd since Cd concentration in 
the leaked brine source is about one order of magnitude lower than 
those of As and Pb. When measured with no -impact thresholds, the 
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Figure 6 | The histograms of the computed plume volumes for pH, TDS, As, Cd and Pb measured with their MCLs and no-impact thresholds at 200 
years. 



SCIENTIFIC REPORT: | 4 : 4006 | DOI: 1 0.1 038/srep04006 



5 



plume volume distributions of As and Pb are similar to those of TDS. 
The volumes of the risk proxies measured with thresholds represent 
the potential areas (or volumes) of contamination (but not over the 
MCLs) in the shallow aquifer. Since the pH plume sizes are the largest 
among the five risk proxies, they will be used as the reference for 
determining the area of review and to design the monitoring network 
for the unconfmed shallow aquifer. 

Implications. An integrated Monte Carlo (MC) simulation under a 
probability framework has been developed in this study for assessing 
risk proxies (plume volumes and probabilities) related to shallow 
groundwater resources in the unconfmed carbonate aquifer. Risk 
proxies are compared to both MCL and no -impact consequence 
thresholds. Differences in risk proxies between the two example 
consequence thresholds (MCL and no-impact) are significant, espe- 
cially for the trace metals. This implies that careful consideration of 
consequence threshold definitions must be made. 

Even though our study is for a hypothetical site with C02 and 
brine leakage, results of our study provide substantial insights into 
the impact of the potential CO2 and brine leakage on shallow 
groundwater quality for the range of parameters studied. Results 
presented here will be used to develop reduced order models for 
shallow aquifers which will be incorporated in lAMs and subse- 
quently used in risk assessment calculations. When more site- 
specific aquifer parameters are obtained from laboratory- scale to 
field-scale tests, calculations of uncertain parameters must be revis- 
ited to refine estimates of all risk proxies. Finally, since the pH plume 
sizes are the largest among the five risk proxies, they will be used as 
the reference for determining the area of review and to design effec- 
tive monitoring network for shallow aquifers. 

Methods 

Previous studies in the Edwards aquifer provide plentiful information about the 
ranges and distributions of the aquifer heterogeneous parameters, such as the aquifer 
thickness, hydraulic gradient, conductivity, and porosity^^"^^. Based on the statistics of 
these parameters, we develop a comprehensive methodology for conducting an 
integrated MC simulation of water and CO2 flow and reactive transport in the 
unconfined aquifer. The integrated MC simulation sequentially uses the uncertainty 
quantification code PSUADE^^, the geostatistical modeling code GEOST^"'^^ (modi- 
fied from GSLIB^'') and the multi-phase reactive transport simulator FEHM^^. 
PSUADE is used to sample all of the uncertain parameters based on the statistical 
analysis results of the field data at or near the study area, conduct global sensitivity 
analysis of the output variables in relation to the input parameters, and evaluate the 
statistical distributions of the pH, TDS, and trace metal plume volumes. GEOST is 
used to analyze the heterogeneous structures of the Edwards aquifer and to simulate/ 
generate the heterogeneous fields for permeability and porosity based on the het- 
erogeneity parameters sampled from PSUADE. 

The multi-phase reactive transport simulator FEHM is used to model CO2 and 
water flow and reactive transport in the unconfined aquifer for each generated het- 
erogeneous model for the unconfined carbonate aquifer. The Brooks and Corey^^'^*^ 
relative permeability is used for water/C02 multiphase flow simulations and the 
related coefficients are adopted from a literature review^^"^*^. The numerical simula- 
tions always start from steady- state flow simulations as the initial conditions and then 
simulate the CO2 multi-phase reactive transport in the aquifer for 200 years. In the 
transport model, it is assumed that the molecular diffusion coefficient is 10"^ m^/s 
and the longitudinal and transverse dispersivities are 800 m and 500 m, respectively. 
For each simulation a post-process is conducted to evaluate the four risk proxies, 
including the plume volumes of pH, TDS, and trace metals, and the CO2 discharge 
rate into the atmosphere at 20 different times, which means that each risk factor will 
have 20 time-dependent values. During the MC simulations we keep checking the 
variance and mean of the risk proxies. After 300 MC runs we find that their variance 
and mean are approaching constants. Then, we stop our MC simulations after 
completing 460 runs when the MC simulations are considered to be converged. By 
taking into account the 20 times as an additional dimension, we obtain the total 
number of outputs of 9200 (realization number (460) times the number of time steps 
(20)). These data are re-organized in a PSUADE input format for global sensitivity 
and risk analysis conducted under a probability framework to evaluate quantitatively 
the chosen risk proxies. 

In order to determine the key flow and transport parameters driving CO2 transport 
behavior in the shallow aquifer, global sensitivity analysis techniques are used for 
investigating input-output sensitivities that are valid over the entire range of the 
parameter sampling. The Monte Carlo simulations provide 460 realizations of input 
transport parameter sets generated with Latin Hypercube sampling and geostatistical 
modeling. Each realization is propagated through the transport simulator to yield the 



output response functions represented by the output variables. Global sensitivity 
analysis entails the comparison of the output distributions to each of the input 
parameter distributions and identifies the most sensitive parameters for each output 
variable. The main effect method^^ is used to quantify the uncertainty and sensitivity 
of the output variables to input parameters. 

The main effect method is a variance-based analysis and it displays first-order 
Sobol' indices for the response surface built from the Monte Carlo simulations^^. The 
essence of this analysis is the statistical measure called variance of conditional 
expectation. The variance-based analysis uses the following equation: 

VCE(XO ='-^t(rj-Yr-^,tti^>J-W, (1) 

;=1 ;=1 i=l 

where, VCE measures the variability in the conditional expected values of output 
variable Y when the input parameter takes on different values, s is the number of 
distinct values of each input parameter and r is the number of replications. N = sr is 
the sample size. The sensitivity of the output variable to the input parameters is 
quantified by Equation ( 1 ) and re-ranked with numbers from 0 to 100 to represent the 
importance of the input parameters'^. 
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